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We develop a systematic approach to the linear-noise approximation for stochastic reaction sys- 
tems with distributed delays. Unlike most existing work our formalism does not rely on a master 
equation as a starting point, instead it is based on a dynamical generating functional describing the 
probability measure over all possible paths of the dynamics. We derive general expressions for the 
Gaussian approximation for a broad class of non-Markovian systems with distributed delay. Exem- 
plars of a model of gene regulation with delayed auto-inhibition and a model of epidemic spread 
with delayed recovery provide convincing evidence of the applicability of our results. 
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Introduction. The theory of discrete Markov pro- 
cesses is well established, and has found applications 
in a variety of disciplines, including biology, chem- 
istry, physics, evolutionary dynamics, finance and the 
social sciences pQ . The standard mathematical treat- 
ment is the chemical master equation, describing the 
time-evolution of the underlying probability distribu- 
tion over states. In simple linear models the master 
equation can be solved exactly [2]. This analytical 
tractability is an exception, although it includes no- 
table examples such as the voter model [3]. The ma- 
jority of Markovian systems can only be analysed us- 
ing approximative schemes, the most naive one is a 
simple mean-field approach. More complex approxi- 
mations are based on the celebrated system-size ex- 
pansion due to van Kampen [4] or on related methods 
such as the Kramers-Moyal expansion [2J [5] . Trun- 
cating these expansions after sub-leading order leads 
to a Gaussian approximation. When linearised about 
a deterministic trajectory this is known as the linear- 
noise approximation (LNA) in chemistry and biology 
4 . Approaches of this type have been applied to a 
wide range of problems [6], and for many model sys- 
tems they reflect the current state of play. Schemes 
going beyond Gaussian order are only currently being 
constructed [?]• 

The purpose of our work is to develop a compre- 
hensive picture of the LNA for interacting-particle 
systems with delay. These are systems in which the 
time-evolution does not only depend on the present 
state, but also on the path the system has taken to 
arrive at this state. Existing work on such models 
include Fokker-Planck equations for delay processes 
8 . Bratsun and co-workers have used approxima- 
tion schemes based on time-separation arguments to 
factorise probability distributions generated by delay 
systems [5] . The system-size expansion to first order 
has been carried out in [5J [TU] in a model with one 
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fixed delay time. Recent work by Lafuerza and Toral 
[TT] has developed these approaches further and con- 
tains the first treatment of systems with distributed 
delays. Most of this existing work is based on the 
formulation of extensions of the master equation for 
delay systems. We here take a different approach 
and choose a generating function description of entire 
paths of the stochastic dynamics [12] ■ This formal- 
ism is originally due to Martin, Siggia, Rose, Janssen 
and Dc Dominicis (MSRJD), and it is not to be con- 
fused with a generating function approach to solving 
master equations. Applying the MSRJD-formalism 
removes the need for a master equation altogether. 
This provides a new perspective on stochastic delay 
systems, and, we think, it allows one to carry out the 
LNA more naturally and systematically. The tech- 
nique we use does not rely on the approximations 
and exclusions required by other methods. As a con- 
sequence we are able to derive a closed and explicit 
Gaussian approximation for a broad class of delay 
models, ready to be applied to problems with delay 
dynamics in a number of fields. 

Generating functional approach to delay systems. 
Consider a reaction system with S types of particles, 
a = 1, . . . , S. The state of the system is characterised 
by n(i) = (ni(t), . . . ,ns(t)), where the integer n a (t) 
indicates the number of particles of type a in the 
system at time t. Assume further that the dynamics 
occurs via R possible types of reactions, i = 1, . . . , R. 
The rate with which each reaction fires is denoted 
by Tj(n). We will consider dynamics in which each 
reaction can result in a change of particle numbers 
at the time the reaction is triggered, and at a later 
time. The latter aspect reflects the delay interaction. 
Specifically we will write Vi yOC for the change in the 
number of particles of type a at the time a reaction 
of type i is triggered. Additionally when a reaction 
of type i fires at time t a delay time r > is drawn 
from a distribution Ki(-). A further change of parti- 
cle numbers then occurs at time t+r, indicated by the 
variables w[ a . This description includes Markovian 
processes, one then has wj a = 0. 
The purpose of expansion methods is to construct 
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Gaussian stochastic differential equations (SDEs) ap- 
proximating the statistics of the reaction dynamics 
[5J These procedures rely on a large parameter, 
N, in most cases a scale setting the number of par- 
ticles in the system. Time is scaled so that reaction 
rates are of order N, Tj(n) = A?Vj(x), and relative 
particle numbers x a = n a /N are introduced. An 
expansion in negative powers of N then leads to an 
effective SDE for x, valid in the limit of large, but 
finite N. Equivalent effective SDEs can be obtained 
using a theorem due to Kurtz [13] . These techniques, 
however, are only applicable for Markovian systems. 

The starting point for our generating function ap- 
proach is a discretised dynamics. Introducing a time 
step A we assume that the number of reactions of 
type i firing at time step t and with a delayed ef- 
fect precisely r € NA time steps later is a Poissonian 
random variable, kj t , with mean Nri[x.(t)]Ki(T)A 2 



We set Ki(r) = for r < 0, and introduce 



(see Appendices A 3 and B 2 ). We will write V(k) for 
their joint distribution, suppressing the dependence 
on x. The generating function for the discrete-time 
process is then given by 

k \ t.a ) 

x JJ<5 [x a .t+A - x a ,t - 0a(k)] ■ (!) 



We have here introduced the source term ip whose 
role is to generate the moments of the {x aj t}- The 
(re-scaled) total change of the number of particles of 
type a at time step t is given by 



6 (k) = N- 1 Y J KtVi 



T>A 



, (2) 



where k i t = X)t>a ^lt- By writing the 8- functions in 
Eq. in their exponential representation, perform- 
ing the average over the {kj t }, keeping only lead- 
ing and sub-leading terms in an expansion in powers 
of A -1 , and subsequently taking the limit A —> 
a continuous-time generating functional is obtained. 
These steps are described in detail in Appendix [B] 
The resulting generating-functional is equivalent to 
the Gaussian dynamics 



= F Q (i,x)+A- 1 /V, 



with (r]a(t)r)p(t')) = B a ^(t,t' , x), and where 



(3) 



ri[x(t)]v i>a + I dTKi(T)ri[yi{t - t)]io£ q 
o 

(4) 



t{t)]Vi t0l 1 



+ dr n[x(t - T)]Ki(T)wl a wlp 



xv itCc w^ t] + ri [ X (t')]Ki(t - t')v it 0W^ a n 



n[x(t)]Ki(t' - t) 
(5) 



Eqs. ([3j |4| |5[) are the main result of our paper, they 
provide general expressions for the Gaussian approxi- 
mation of a wide class of delay systems [HI [T5] . This 
result is slightly stronger than the LNA [2] , the latter 
can be obtained from Eqs. (|3]|4j[5]) by a straightfor- 
ward linearisation (see Appendices A 6 and B 4 1 . In 



the following we will demonstrate the applicability of 
this approach. We will use our results to compute the 
spectra of noise-induced quasi-cycles [H] in a model 
of gene regulation and in a model of epidemic spread, 
both with delay interactions. We stress that the char- 
acterisation of quasi-cycles is only one of many ap- 
plications of the Gaussian approximation or LNA of 
delay systems. 

Application to a model of gene regulation. Delay re- 
actions play an important role in gene regulatory pro- 
cesses. Delays in transcription and translation pro- 
cesses are considered a potential mechanism for os- 
cillatory behaviour in somitogenesis, these cycles are 
thought to give rise to spatially heterogeneous cellu- 
lar structures [TTJ [TS] . Models of these processes have 
traditionally focused on differential equations, see e.g. 
[18]. It is only more recently that intrinsic noise has 
been included in models of gene regulation [51 H"9"] . 
This is based on the observation that particle num- 
bers in gene regulatory systems can be small, making 
deterministic approximations inadequate 20J . For 
example noise-driven quasi-cycles go undetected in 
deterministic models |16| . Existing theoretical analy- 
ses are limited to models with constant delay periods 
[10] , we note recent advances [TTJ . Our result for 
systems with distributed delay provides a systematic 
theoretical framework, and we apply it to the simple 
model of gene regulation described in 18,19 . We will 
consider two types of particles, mRNA molecules, de- 
noted by M, and protein molecules, P. The stochas- 
tic dynamics are given by 



M 
P 
M 



g(n P ),K(r) 



M + P, 
M. 



(6) 



The first two interactions correspond to degrada- 
tion of mRNA and protein, respectively, the con- 
, stant model parameters fiM and /zp describe their 
degradation rates. The third interaction describes 
the translation of mRNA into protein. Finally, the 
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FIG. 1: (Colour on-line) Power spectra of quasi-cycles 
in the gene regulatory model [181 119] with uniformly dis- 
tributed delays over the interval [18.7 — ft/2, 18.7 + ft/2] 
minutes. Lines are theoretical predictions within the 
LNA, markers from simulations using a modified next- 
reaction method and averaged over 700 realisations. 
Parameters are om = ap = 1, /ip = /ijvj = 0.03 (all with 
units min" 1 ), P = 10, h = 4.1, N = 500. 



fourth interaction represents the transcription pro- 
cess, within the model effectively the production of 
mRNA. This process is suppressed by the presence of 
protein molecules, as reflected by the Hill function, 

g(n P ) = a M [l + [n P /{P Q N)] h ]~\ where h and P 
are constants. The double arrow indicates a delay 
reaction. In this particular model the reaction has 
no effect on particle numbers at the time t it is trig- 
gered, but only at a later time t + r, where r is a 
distributed delay time drawn from K(t). The reac- 
tion rate depends on the number of proteins at the 
earlier time, np(t). Earlier works [lOj [19] focus on 
the case in which K(-) is a ^-distribution, and ex- 
clude distributed delays. Applying our general result 
above (see Appendix [C| for details) we find 



x M (t) 
i P (t) 



dr K( T )f[x P (t-T)} 

-HttXM{t)+N~ 1,2 VM$), 
a P x M (t) - fipx P (t) + N-^ 2 VP {t), (7) 

where f[x P (t)} = [1 + (x P (t)/P ) h }-\ and 



(VM(t)r, M (t r )) = 



(VM(t)vp(t r )) = 



a M dr K(r)f[x P (t-T)} 



6(t-t'), 



+HMX M (t) 



[a P x M (t) + jj,pxp(t)\ 6(t - <'), 
0. (8) 



The Gaussian noise components r]M,iJP have no cor- 
relations in time, as expected for a dynamics in which 
each reaction changes particle numbers only at one 
single time. A more complex case will be studied be- 
low. In the deterministic limit, N oo, Eqs. Q, 



are typically found to have a fixed point, (x* M ,x* P ). 
A systematic expansion, xm = x* M + A^ 1 / 2 ^/, and 
similar for xp, then leads to the LNA: a pair of 
linear Langevin equations for the fluctuation vari- 
ables £a/ and £p. A straightforward calculation fol- 
lowing the lines of [TB] then allows one to compute 
the power spectra of noise-induced cycles, Pm{u) — 

{\£,m{u )\ 2 ^J , and similarly for the protein (see Ap- 
pendix |C 2\ . Results for a uniform distribution of 
delay times are shown in Fig. [T] and are confirmed 
convincingly in numerical simulations. 
Application to a model of epidemic spread with 
delayed recovery. We consider a variant of 
the susceptible-infective-recovered (SIR) model with 
birth and death [55] . Specifically the model describes 
a population of a fixed size of individuals, each of 
which can be in one of three states, susceptible (S), 
infective (I) or recovered (R) . Infection occurs via the 

process S + 1 — > 21, and the newly infected individ- 
ual may recover (/ — > R) at a later time, where the 
delay is drawn from a distribution H(-). However, 
all individuals are subject to a birth-death process, 
occurring with rate fx, and in which an individual 
dies and is immediately replaced by an individual of 
type 5*. This is a commonly used simplification, en- 
suring a constant population size "22 . This setup 
implies that a newly infected individual may die and 
be replaced by an individual of type S before its des- 
ignated recovery time is reached. This is illustrated 
in Fig. [2] Assume an infection occurs at time t. 
One may think of the subsequent dynamics as fol- 
lows: at the time of infection, a designated time-to- 
recovery, t, is drawn from H(-). At the same time 
a designated time-to-removal, s, is drawn from an 
exponential distribution, E(s) = /ie _Als . There are 
then two possible subsequent courses of events: (i) 
If t < s recovery occurs before death, the recovery 
process completes at time t + r, and the infective in- 
dividual is replaced by an individual of type R. The 
death event is discarded. The probability for case (i) 
to occur is x = dTH(r) ds E(s). Conditioned 
on this sequence of events, i.e. if recovery occurring 
before death is a given, the time-to-recovery follows 
the distribution K(t) — x^-HXt) ds E(s). Case 
(ii) describes the opposite situation, s < t, occurring 
with probability 1 — x- In this case the newly infected 
individual dies before the designated time of recovery, 
and we have a reaction of type / — > S at time t + s. 
The conditional time-to-removal, given that case (ii) 
is realised, is Q(s) = (1 - x)~ 1 £^(s) dr H{t). We 
can summarise the reaction scheme as follows 



1? 

S + I 



/-' > 

X/3 



S, 

2I;im R , 



s+i {i =n fi 2Lim S . 



(9) 



The notation for the second reaction channel, occur- 
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b) Individual recovers 

S + I I 
i ♦ 

21 R Discarded 



H(r) 
E(s) 



~1~ 

t + T 



t + S 



c) Death occurs before recovery 



5 + / 
* 

21 



Discarded 



t + s 



t + T 



FIG. 2: (Colour on-line) Possible sequences of events 
when a reaction with delayed recovery is triggered. A 
time-to-death, s, and a time-to-recovery, r, are drawn 
from the appropriate distributions (panel a)). Depending 
on the outcome recovery or death may occur (panels b) 
and c) respectively), the remaining event is discarded. 



ring with rate 72(n) = fix n S n i /N , indicates that 
one particle of type S is converted into an / at 
the time the reaction is triggered, and that an in- 
dividual of type I is converted to R at a later time 
t + T, where r is drawn from the distribution K (•). 
Similarly, the third reaction channel fires with rate 
T 3 (n) = /3(1 — x) n s n i /N , and results in an event 
S + 1 -+ 21 at the time the reaction is triggered, and 
then in an event of type / — > S at a later time t + s, 
where s is drawn from the distribution Q(-). 
Applying the general result above we find (with S = 
n s /N,I = m/N), 



S(t) = -pS(t)I(t) + n{l-S(jk)-I(t)) 

+0(1 - X ) I dt' Q(t - t')S(t')I(t') + N-WrisM, 



I(t) = 0S(t)I(t) - (3 / dt'S(t')I(t') 




xK(t-t') + {\-x)Q(t~t') 



N- 1 / 2 ^). (10) 



Unlike in the above model of gene regulation, the 
noise is now correlated in time. Expressions for the 
correlation matrix arc lengthy and reported in Ap- 
pendix [d] 

Recent theoretical work has studied SIR-models in 
which individuals progress through a series of L in- 
fectious 'stages', 1% -+ I2 -+ ■ ■ ■ -+ II — > R, at rate 
jL before they recover (or die along the way) [35] . In 
our formalism this is equivalent to a model in which 

H(-) is a T-distribution, H(t) = i j^yT L - 1 e-~< LT . 
To make contact with the results of [55] we use Eqs. 
( 1 1 to compute the power spectra of noise-driven 



quasi-cycles about the fixed point of the deterministic 
limiting dynamics (see Appendix D 2 1 . Results from 
the theory and from simulations are shown in Fig. [3] 
We have carried out simulations of the staged model 



FIG. 3: (Colour on-line) Power spectra Pi{oj) for the SIR 
model with delayed recovery. Lines show results from 
the LNA for the staged model (SM), see [22j |23], and 
for the delay model (DM). Markers are from simulations, 
averaged over 800 independent runs. Model parameters 
are /3 = 10.56, fi = 4.81 ■ 10 _3 ,7 = 1. System size is 
N = 10 6 . Inset: Results for L = 4 and ju = 4.81 ■ 10~ 2 . 



(SM) and of the delay model (DM). Both give identi- 
cal results, confirming the equivalence of the SM and 
the DM with a T-distributed kernel. Simulations of 
the staged model are less costly than those of the de- 
lay model. The analytical calculation of the results 
in Fig. [3] is more demanding in the staged model, 
as it involves a larger number of particle types. It 
should be noted that the staged model is limited to 
T-distributed recovery times, whereas our approach 
is more general. For sufficiently small values of the 
death rate, /z, there is no noticeable difference be- 
tween the predictions of our approach and the result 
of [55] (see main panel). The result of [55J is based 
on an expansion in /j, and deviates from simulations 
when the small-^t approximation is not justified. Our 
theory does not rely on such approximations, and de- 
scribes simulation results accurately in such cases (see 
inset of Fig. [31). 

Conclusions. We have presented a comprehensive ap- 
proach to the LNA for stochastic dynamics with dis- 
tributed delays. Our calculation is based on a gener- 
ating functional, rather than a master equation. In- 
stead we focus on probabilities to observe entire paths 
of the dynamics. This makes the approach suitable 
for non-Markovian systems, and we are able to de- 
rive general expressions for the Gaussian approxima- 
tion of a broad class of processes with distributed 
delays. The applicability of our results is demon- 
strated through the computation of power spectra of 
noise-driven cycles in delay models of gene regulation 
and of epidemic spread. We expect that the gen- 
eral expressions we have derived will be of use for a 
variety of phenomena in the biological and physical 
sciences, and indeed in other areas where individual- 
based models with delayed interactions are relevant. 
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Appendix A: System-size expansion without a master equation: Markovian systems 

1. General remarks 

In this section we describe an alternative method with which to obtain results of the well-established expansion 
methods for Markovian systems. These expansions are usually carried out starting from a master equation 
describing the underlying stochastic process, and then following the procedure originally first proposed by 
van Kampen [3], or alternatively the steps of the Kramers-Moyal expansion [2J. We here choose a different 
approach, and use a moment-generating function as the starting point. This technique is originally due to 
Martin, Siggia, Rose, Janssen and De Dominicis (MSRJD) [12]. 

The central object underpinning the generating function approach is the probability measure of possible 
dynamical paths of the system, and not the probability of finding the system in a given state at a specific 
time. We will first discretise time in steps of duration A, and keep A finite while we perform the average over 
paths and while we carry out the system-size expansion at the level of the generating function. The limits of 
large N and A — » are therefore decoupled, see also [23] for a similar approach in a different context. It is 
only at the end that we take the the time step to zero, A — > 0. This then results in a generating functional 
for a Gaussian process in continuous time, equivalent to the Gaussian dynamics one would have obtained 
from a direct Kramers-Moyal expansion of the underlying master equation or from Kurtz' theorem [13:. This 
multiplicative Gaussian process can then be linearised straightforwardly. For Markovian processes our method 
reproduces the known results of the linear-noise approximation. 

The generating function approach we take here is related to, but different from the path-integral formalism 
introduced by Doi and Peliti |25) . The latter starts from a master equation, which is not required in the 
MSRJD-formalism. For Markovian systems the two approaches can be shown to lead to the same results, see 
for example [25] for a discussion. We also note that the outcome of van Kampen's system-size expansion for 
Markovian systems can be obtained in the Doi-Peliti formalism [27] . 

The true strength of the generating function method becomes more transparent when we turn to non-Markovian 
delay dynamics. In such cases there is no closed master equation, and so we feel the generating function 
approach, which focuses on entire paths of the dynamics, is the most transparent technique to deal with delay 
dynamics. 

2. A Markovian example - model definitions 

We first illustrate the method using a specific example of a Markovian system, defined by the following set of 
reactions 

-±> A, 
A + B 3B, 

B A 0. (Al) 

This notation indicates that the three reactions occur with rates 

TLITi 

T 1= N, T 2 = a—, T 3 = bm. (A2) 

We have here written n for the number of particles of type A in the system, and m for the number of B- 
particles. The quantity N is a model parameter, and sets the scale of the system size. In a continuous-time 
setting reactions occur as exponential processes, the quantity Ti[n(t) , m(t)]dt represents the probability for a 
reaction of type i to occur in the time interval (t, t + dt). 

3. Discretisation of time and the generating function 

In order to define the generating function we will discretise time into time steps of duration A. This requires 
further qualification of the above reaction rates. We will be looking at a discrete time process, with paths 
{n t ,m t } = . . . , [nt-A,mt-A), (n tl m t ), {n t+Al m t+A ), .... 
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The dynamics of this system is determined by how frequently each of the three reactions above fire. We 
will write fc^t for the number of reactions of type i — 1,2,3, which occur between time t and time t + A. 
Consistent with the above reaction rates we will assume that k it is a Poissonian random variable with mean 
Ti(nt,m t )A.. This amounts to a discretisation of the above continuous-time exponential process, similar to 
what is considered in the context of the algorithm with so-called r- leaping j^Sj. The starting point of our 
calculation is the Martin-Siggia-Rose-Janssen-De Dominicis (MSRJD) moment-generating function 1 [12], see 
also I2T 



X t+ A - X t 



x exp iA ^2 Vl^tXt + <ftVt 




2k 



yt+A - yt 



2,1 



N 



N 



(A3) 



We have here introduced the variables x t = and y t = and we have written flx = Y[ t dx t , and similar 
for Dy. We discuss the absence of any Jacobian determinants in the generating function below in the context 
of the more general model with delay (see Sec. [B]). It should be noted that t comes in integer multiples of 
the time step A, i.e. we have t = IA, where £ = . . . , —2, —1, 0, 1, 2, ... . Our notation always implies this 
convention, it is important to keep this in mind when it comes to objects such as iptX t , a short-hand for 

The quantity 'P(k) describes the joint probability distribution of the {fci,*}, we note that the statistics of 
k it depends on x t and y t - This dependence has been suppressed in our notation. The source terms xp = 
{tpt} and tp — {(fit} finally have been introduced as per normal procedure [T^]. Taking derivatives with 
respect to these source terms generates moments of the variables {xt} and {yt}, for example one has (xt) — 

ip=cp=o 



where 



denotes the average over all possible paths of the system, 



</[x,y]>=^ f DxDyV(k) 

1, 



IP 



Xt+A ~ Xt 



k\, t 
N 



k2,t 
N 



yt+A - yt 



2k 2 ,t 
N 



k 3 ,t 
N 



/[x,y]. (A4) 



4. Average over trajectories and system-size expansion 



Starting from Eq. ( A3 ) we first write the delta- functions in the integrand in their exponential representation 
and obtain 



k J 

x exp li 



n 



dx t dx t dy t dy t 
2tt 2tt 



V(k) exp iA ^ {ij; t x t + tp t y t } 



kit , k 2 . t 

X t+ A-X t - W + W 



Vi 



2k 



yt+A - yt 



2,t , K 3,t 



N 



N 



Next we collect terms and perform the average over the {fcj,t}. The relevant terms are 

^■P(k)exp jr;^2{- k i,txt + k 2>t (x t - 2y t ) + fc 3 ,t& J . 



(A5) 



(A6) 



The distribution over the {k i t } factorises, 'P(k) = ]X Tit Pi,t{ki, t)- The individual factors, P i t (-), are here 
Poissonian distributions with parameters A^t = Ti(n t ,mt)A. 

As a specimen and to demonstrate the procedure, we now carry out the averaging over one of the fc^t, 
specifically we choose k 2 < for a fixed time step t. We have 



k 2 ,t 



k 2 J 



N 



k 2 ,t 



\X-2_te 1 n 



-in 



fc2,t 



k 2 , t \ 



exp A 2 ,t + A 2 , t e i 



(A7) 



1 We will use the term 'generating function' for systems in discrete time, and 'generating functional' when time is continuous. 



7 



The next step is to expand this expression in powers of N 1 . We obtain 



exp (-A 2 , + X 2tt e^) = exp ~ ^J 1 ^ 1 + ^ * 0(N' 3 ) 



(A8) 



Keeping only terms of leading and sub-leading orders we therefore have 



A 

Ai,t M — A 2 ,t 2,t —A3 ,t 3,t 



l fe 3,t 



exp ( — y^(-fei,tJt + fc 2 ,i(£t - 2 2/t) + fc 3. 



exp 



i~(-Ai,t + A 2 ,t) + i|(-2A 2 , 4 + A 3 , t ) - ^(A M + A 2)t ) - -f, (4A 2 ,* + A 3 , t ) 



TV 2 



2A 2 



(A9) 



and so we arrive at the following expression for the generating function 



n 



dxtdxt dy t dy t 
2tt 2tt 



exp i^x t (x t +i - x t ) + Vt(vt+i ~ Vt) exp iA^x t (-ri lt + r 2 ,t) 



exp ^* A E yt(- 2r 2,t + r 3 ,t)^ exp E ( r M + r2 -*) ~ Vt{^ r 2,t + r 3 ,t) - 4x f y t r 2 , t }^ 

x exp ^iA ^ [V'txt + . (A 



We have here introduced the quantities r^t- These are shorthands for ri(xt, yt), which in turn are defined via 
written \ i>t = Ti(n t ,m t )A = NAri(x t ,yt), keeping in mind that the rates Tj are of order N. The quantities 
fi,t — r i(xt,yt) are therefore of order N°. 



Expression (A10) can be re- written as 



z[1>M 



n 



dx t dx t dy t dy t 
2tt 27r~~ 



exp ia > xt 



X t+ l - x t 
A 



ri.t + r 2 ,t 



exp ^iA ^ 



Ih 



Vt+i - yt 
A 



2r 2 , t + r 3 , t n exp UA^[V>, 



t a; t + iptyt] 



exp i -^7 ^ {x~tBx,x,t,t'Xt' + ytBy,y,t,t'Vt' + %XtBx,y,t,V Vf } 



(All) 



with 
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®y,y,t,t< 



(4r 2>4 + r 3 1 ) 



A ' 
A ' 



= -2r 2 



' A ' 



(A12) 
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5. Continuous-time limit 

Taking the limit A — > we find 

Z[i>M = J Vx.V9.VyVyexp(i J dt x(t)[x(t) - n[»(t),y(t)] + r 2 [x(t),y(t)}} 
exp (i J dt y(t) [y(t) - 2r 2 [x(t), y(t)} + r 3 [x(t), y(t)}]\ 

xexp ("2^ / dtdt ' {x( t ) B xA t i tl )x( t ') + y( t ) B v^ tr W) + 

x exp U I dt \i>{t)x{t) + ip(t)y(t)} ) , (A13) 



where the (path-) integral J 2?x. . . runs over all continuous-time paths {x(t)}, and similarly for J T>y . . . and 
the auxiliary variables [12] ■ We have also introduced 

B XlX {t,t') = r 1 [x{t),y{t)]+r 2 [x{t),y(t)]S{t-t'), 
B yiV {t,t') = 4r 2 [x{t),y(t)]+r 3 [x(t),y(t)]6{t-t'), 

B XlX (t,t') = -2r 2 [x(t),y(t)]S(t-t'). (A14) 

We have here used the correspondence A~ 1 5 t .t' *H> S(t — t') between the Kronecker-J for discrete arguments, 
and the Dirac (5-function for continuous arguments. This correspondence can easily be verified using the 
correspondence AJ2t ft ^ J dtf(t), as well as / dt'S(t — t')f(t') = /(i), and the observation that J2t' $t v ft' — 
A&(A"V)/. = / ( . 



The expression in Eq. (A13) is recognised as the generating functional of the following dynamics 



x(t) = r x [x(t),y{t)]-r 2 {x{t),y{t)} + ~ V {t), 



y(t) = 2r 2 [x(t),y(t)]-r 3 [x(t),y(t)] + ^={(t), (A15) 



with white Gaussian noise variables rj(t) and (both of mean zero), and with correlations 

( V (tW)) = {ri[x(t),y(t)] + r2[x(t)Mm^-^ 
(C(«)C(i')> = {ir 2 [x(t),y(t)}+r 3 [x(t),y(t)}}S(t-t% 

(V(t)m) = -2r 2 [x(t),y(t)]S(t-t'). (AW) 

This is the result one would have obtained from a direct Kramers-Moyal expansion of the master equation of 
the system [2] or by applying Kurtz' theorem [13] . For later convenience we write Fi(x, y) — r\(x, y) — r 2 (x, y) 
and F 2 (x,y) = 2r 2 (x,y) — r 3 (x,y). In our specific example one has r\[x(t), y(t)] = 1, r 2 [x(t), y(t)] = ax(t)y(t) 
and r 3 [x(t),y(t)] = y(t). 

6. Linear-noise approximation 



The result obtained in the previous section, Eqs. ( A15[ A16), describes a process with multiplicative noise, 



and is a slightly stronger result than the linear-noise approximation. The latter is obtained by writing 
x (t) = x°°(t) + N- 1 / 2 ^) and y(t) = y°°(t) + N' 1 / 2 ^), where the deterministic trajectory x°° (t) , y°° (t) 
is the solution of i°° = Fi (x°° , y°° ) , y°° — F 2 (x°°, y°°). The superscript 'oo' indicates that this deterministic 
trajectory is obtain ed in the thermodynamic limit, N — > oo. 



Inserting into Eqs. ( A15), and expanding in powers of N 1 I 2 , one then finds 



■ dF^x^T) M dFiix^T) 1/2 

& = ^ + ^ycc ^ + N ^ 

fa- dx^ & + Qy^ &+ N 1*, (A17) 
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where 

(li(t)li(f')) - {n[x°°(t),y°°(t)}+r 2 [x°°(t),y°°(t)}}5(t-t'), 

{m{t)m{t')) = {4r 2 [x°°(t),y°°(t)}+r 3 [x°°(t),y°°(t)}}6(t-i;), 

{ m {t)m{t')) = -2r 2 {x°°(t),y°°(t)]S(t-t>). (A18) 
The noise is now linear (additive). 



7. General theory for Markovian processes 



We now consider a general Markovian model with S types of particles, a = l,...,S and M reactions, i = 
1, . . . , M, again in discrete time. We will write n Qjt for number of particles of type a in the system at time 
t, and we introduce x = (xi, . . . ,xs), where x a — n a /N. We assume that the reactions occur with rates 
Ti = ANr i} where the Ti is of order TV , and may depend on the state of the system, ri = r^x). We will 
assume that the occurrence of a reaction of type i will result in a change of the number of particles of type a 
by i.e. the number of particles of type a will be n a j+A — n a , t + h^Vi^ if k itt reactions of type i occur 
between time t and time t + A. 
The generating function is now given by 

Z W\ = X] / ^ x ^( k ) n 5 l x a ,t+A - x a ,t - kt,t ^' a ) cxp I ip a , t x a ,t ) 

k lt,a \ i ) J \ t,a J 



1, J 



n 

,a,t 



dx f_ dx ^ 



V(k) exp I i X x Q , t I x Q , t+A - x att - ^ 



N 



exp iA X V 



The terms containing a given Poissonian variable k itt are now 



fci,, 



exp ^-^\i,ty^x a ,tVi,<x ~ ^?.^i,t (^2x<*,tVi,c}j 

exp j -JjKt ^2%a,tVi, a - ^pKt ^2xa,tVi t aVi,pXpj + . . . J . 



This leads to 



n 



dxex^dxtx^ 
2^ 



A 



exp iA X 5? a ^ t 



%a,t+A — X a .t 



x exp ( -— X X x a ,tB a ,i3,t,t<xi3,t> + ■■■ ] exp ( zA X tp a ,tx a ,t ) , 



i,t' q,/3 



up to and including terms of sub- leading order in TV 1 . We have here introduced 



In the limit A — >• we find 



(A19) 



(A20) 



(A21) 



(A22) 



Z[V>] = J PxPxcxp y dt x a (t) ^x a (t) -^r,[x(()]« l!a jj 

^ J dtdt'^x a {t)B at p(t,t')xp{t') + ...J exp^X y Va(*)a:a(t)^ , (A23) 
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where 

Ba,p(t, t') 



ri[x(t)]vi, a Vi,^J 5(t - t'). (A24) 



This in turn is the generating function description of the dynamics 

x a (t) =Y,n[x{t)Ka + N- 1 / 2 ri a (t), (A25) 

i 

where the {i]a{t)} describe Gaussian white noise of mean zero and with correlations 

(V a (t)vp(t')) = fe«i,a»i >j sr<[x(t)]^ 5(t-t') (A26) 

across particle types. Again this is the result one would have obtained from a Kramers-Moyal expansion [2J, 
or from a direct application of Kurtz' theorem |13j . 

The linear noise approximation is obtained following the steps outlined above. One writes x(i) = x°°(t) + 
N' 1 / 2 ^), and finds 

Ut) = £ J a Armpit) + Tla(t) (A27) 

P 

where J a ,g = dF a /dxp, and where the replacement x(i) — > x°°(i) is implied in the noise correlator given in 



Eq. (A26). This makes the noise additive. We have here written F a (x) = ^\ ri(x)v i a . Eq. (A27l represents 
the result one would have obtained from a direct application of van Kampen's system-size expansion [3]. 



Appendix B: General theory for delay systems 



1. Model definitions 



We will next look at a discrete-time system with delay interactions. We will again assume that there are S 
types of particles, a — 1, ... ,S, with particle numbers n a j, and that there are M reactions, i = 1, . . . ,M. 
As before we write x a = n a /N. In the continuous-time model reaction i is taken to fire with rate Tj — Nri, 
where r.; = r.i(x) ~ O(N ). The immediate change in particle numbers at this time will be denoted by Vi^ a as 
before. 

A delay reaction firing at time t can affect particle numbers at later time steps. Throughout this paper we 
assume that any instance of a reaction will only induce changes of particle numbers at at most one single 
later time (this is not a severe restriction and will be the case for most applications). This time is a random 
variable, drawn from an underlying probability distribution. We will further restrict the discussion to models 
in which the distribution of delay times, r, of a given reaction does not depend on the state of the system 
when the reaction is triggered. Again this is not a severe constraint, and is usually fulfilled (generalisation of 
our theory to models in which this is not the case are possible) . To be more specific we will assume that once a 
reaction of type i is triggered at time t an immediate change of particle numbers by Ui jCt occurs (a = 1, . . . , S). 
Additionally a delay time r > is drawn from a distribution Ki(-), and a further change of particle numbers 
by wj a then occurs at time t + r. The appropriate normalisation of K(-) is J °° drK^T) = 1. The rate with 
which a reaction of type i with a delay time of precisely r fires at time t is hence Tj [n^Jifj (t) . The probability 
to see a reaction of type i fire during the time interval from t to t + dt and with a delay time in the interval 
from t to t + dr is Ti(n(<)) Ki(r)dtdT . 

This change in particle numbers at the later time can in principle depend on the drawn delay time r, hence the 
superscript in wj aM , even though this is not usually the case in most applications. We allow for this possibility 
when we develop the theory, in the two applications we study in the main paper the wj a are independent of 

T. 

We also point out that our formalism allows one to include reactions without delayed effect. These will simply 
have wj a = 0. In absence of any reactions with delay (wl a = for all i, a, r) we recover the Markovian case 
discussed in the previous section. 
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2. Generating function approach 

In the discrete-time setting (with time step A) the number of reactions of type i triggered (with any delay 
time t) at time step t will be a Poissonian random variable ki t t with mean TjA = NAri(x t ). The number 
of reactions of type i with a delayed effect t = £A time steps later will be a Poissonian random variable k T it 
with mean (ATi) x (AKi(r)). The second factor, A x Ki(r), is the probability that a delay time in the time 
interval (r,r + A) is drawn in the continuum model. We have kit = kj t . For later convenience we will 
introduce pi jT ~ A x Ki(r), so that kj t becomes a Poissonian random variable with mean A[ t = NAri(^c t )pi,T- 
The variable k^t follows a Poissonian distribution with parameter \ i t = ^ T XJ t = _/VArj(x t ). 

The MSRJD generating function for this process is given by 2 



z[i>] = E / ^ x w 



II 5 [x^t+A-Xoct-^Kt^ -ee 



i t>A 



k lt-r W i,< 

N 



exp liA^ip 



(Bi) 



In order to keep the notation sufficiently compact we will simply write ^ T ... in the following instead of 
S r >A • • • ■ The constraint r > A is always implied. 
Given that ki t t = X) r >A ^It this can be written as 



zty] = E / D *- V M 



X a ,t+A ~ Xa,t 



N 



E 



kJ,t-T W i,a 



E 



n 



dx q _ ^ dx q 
2^T~ 



7>(k) exp i ^ ®a,t %a,t+A - X a .t ~ ^ 



exp ^A^ a , t i a ij . 



exp iA^tp a . t x ay 



k\ t Vi^ a + kJ t _ T wJ c 
N 



(B2) 



The terms containing the Poissonian variables k are given by 



ex P [ -Jj E Kt E {%»,tVi, a + ^,Hr<«} 



i.t.r a 



eXP ( ~N E E \ KtVt.a + E A « r ,t-r< C 



exp 



LeeI 2 ^ 



2iV 2 



x exp 



iv2 E E i E 



27V 2 



t a,j9 



\,tVi, a Vi,/3 + ^2 X lt-T w i, a w i,f: 



X0,t 



(B3) 



We have here written (. . .) k for the average over the {kl t }. Putting the pieces together, and replacing 



2 We note that potential functional determinants in this expression are field-independent and can hence be discarded, see also 
26 29-31 . If we introduce the short-hand X at+ & for the expression inside the delta-function in Eq. IBlL then the relevant 
Jacobian is given by ^T a n t,t' = dX a ,t/dxp t i. Keeping in mind that the statistics of kj t only depend on the state xt, but 
not on states, x t / at times t' > t, one then notes that J a stt' = f° r > *i due to causality. One also has 3 a ,p,t,t = &a,p- 
The matrix J is triangular with respect to the indices t and t' , with unit diagonal elements. This makes its determinant 
field-independent and equal to unity. 



NAn(x t ), and XJ t = N Ari(x. t ) p i<T , we have 



m = j 



I] Xa '^ a ' t exp | i ^ %a,t {x a ,t+A - X at t) J CXp | lA ^ Va.t^a.t 
J \ a,i / \ t.a 

ri(xt)vi, a v it p + ^ rj(xt-T)Pi,Twl a wl 

T 

^2^2{ 2x ~a,t ^2n(xt)Pi,TVi, a wlp Xfi.t+T 



A 

2N 



t a,/; 



x exp 



A 



This can be written as 



zm = J 



n 

.ct,t 



2tt 



exp iA ^2 

Xa,t 



X a ,t+A ~ Xa,t 

A 



- f a ,t ) exp I iA ^ TP 



t,a 



exp -—^2^2x a jB a ^j,t>xp,t 



t,t> a,P 



where 



fa,t = X | r i( X t) w i,a + ^2 r i(x t - T ) Pi^wJ^ , 



B, 



a,/3,t,t< 
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N I A 



ri{x t )vi, a Vi,i3 +^2r i [x t - T }p i , T wl a wl l3 



, ST^ [ / sPi,t'-t (t'-t) , / ,Pi.t-t' (t-t') 



3. Continuous-time limit 



Taking the limit A — > we find 



Z[tp] = j 2?x2?xcxp ^A^ J dt x a (x) [x a {t) — F a (t,x)]j exp J ^ V'a(*)a ; a(i) 

xexp I --j- ( dt dt' ^2x a (t)B a ^(t,t' ,x)x p (t r ) J , 



where 



F a (t,x) = |rj[xft)]^, a + J dr n[x 



(t-T)]Ki{T)wT\, 



B aJj (t,t',x) = h(t-t')J2 



ri[x(t)}vi ta Vi^ + / dr r;[x(i - T^Ki^wJ wJg 
o 



+ ^2 [rMt)]Ki(tf - t)vi, a w? fi t] + ri[xtf)]Ki(t - i'K/5»\ 
This corresponds to the dynamics 

x a =F a {t,x)+N- 1 ' 2 T la , 



(t-t') 
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where 



ri[x.(t)]vi^ a Vi,f3 + dr r,[x(f - T^Ki^wJ^wJ^ 



+ J2 [n[x(t)]Ki{t' - t)vi, a w? p t] + r t [x(t')]K t (t - t')v^w^ a °] . (BIO) 



4. Linear-noise approximation 

One proceeds again by decomposing x(t) = x°°(t) + N^ 1 / 2 ^), where x°° is the deterministic trajectory, i.e. 
the solution of 



Substituting in Eq. ( B9 1 we find 



x a — F a (i, x 



k(t)= / dt' J(t-t', X °°)Z(t)+r,(t), 



(Bll) 



(B12) 



where we again imply the substitution x — > x°° in the correlator of the noise, see Eq. (BIO I. The Jacobian is 
defined by 



J a)/3 (i,t',x) = 



(B13) 



and depends only on time differences, t — t' . 



Appendix C: Model of gene regulation with delay interaction 
1. Gaussian approximation 

The delay model of gene regulation discussed in the main paper (see also [HI HH]) describes two types of 
particles, a = P, M. We will write 



x P 

X M 



(CI) 



In the model there are four reactions, i = 1, . . . , 4, one of which is a delay reaction. The reaction rates are 

Vmxm 



r[x(i)l = ^ pXp 

ot M f[xp{t)] 



(C2) 



with f[x P (t)} = [1 + {x P {t)/P Q ) h ]- x . 

The effect of the reactions at the time they are triggered are described by the stoichiometric coefficients 




(C3) 



where each row represents one of the reactions, and where the first column is the change in number of protein- 
molecules and the second column is the change in the number of mRNA molecules. 

There is only one delay reaction in this model, i = 4. The corresponding distribution of delay times is 
Ka(t) = K[t), and the delayed effect on particle numbers is given by w^m = 1. All other uii. a vanish. 
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Applying the general result discussed in the main paper [Eqs. (3, 4 , 5)] we find 

F P (t,x) = a P x M (t) - fJ, P xp(t), 

F M (t,x) = a M f dt'K(t-l/)f[xp(t')}-HMXM(t), 

J — OO 

for the drift terms, as well as 



(C4) 



Bp, P (t,lf,x) = (a P x M (t) + ^px P {t))5(t-t r ), 
B PtM (t,t',x) = fl Jtf|J >(t,i / ,x) = 0, 

%m(M'x) = (»MX M (t) + a M [ dt"K{t-t")f[xp{t")])5{t-t'), (C5) 



for the correlator of the noise rj. 

2. Linear-noise approximation and spectra of noise-driven quasi-cycles 

For the parameters discussed in the main paper the deterministic dynamics, obtained in the limit N — > oo, is 
seen to converge to a fixed point (x* P , x* M ). The stochastic system exhibits noise-driven quasi-cycles. In order 
to predict the spectral properties of these cycles we follow the standard procedure outlined above, see also 
[IB] , and write xp — x* P + £,p/\fN and xm = x* M + £m/VN. Substituting this in Eq. (7) of the main paper, 
we then perform an expansion in powers of N~ 1/>2 . One obtains 

€(*)=/ dt^t-f^tfO + irt*). (C6) 

J — oo 

where J(t — t',x*) is the (functional) Jacobian of the deterministic dynamics, evaluated at the deterministic 
fixed point. Specifically we have J a ,/3(t, t' , x*) = S f°^^ 



, i.e 

FP 



Jp,p{t-t',x*) = -np5{t-t'), 

J P , M (t-t',x*) = a P 8{t-t'), 

J MtP (t-t',x*) = a M K(t-t')f'(x P ), 

JM,M(t-t',x*) = -fi M 8(t-t')- (C7) 

with f'(xp) the derivative of f(xp) with respect to ip. We note that these matrix elements are time-translation 
invariant, i.e. they are functions of t — t' only. We also point out that K(t — t') — for t' > t. 



Carrying out a Fourier transform (with respect to t — t 1 ) of Eq. (C6) we have 

MH€(w) = rj(u), (C8) 

where 

M(w) =iu)L-l{u,x*), (C9) 
and where / is the 2x2 identity matrix. For this model we have that, 

M(lj) = ( , , V (CIO) 

At the deterministic fixed point we have 

Bp P (t-t', x *) = (a P x* M + fipxp)5(t-t'), 
B P>M (t-t'X) = B MiP (t,l/,x*)=0, 

B MM (t-t',x*) = (fXMX* M +a M f[xp)})S{t-t'). (Cll) 
We will denote the Fourier transforms of these matrix elements by B a ^{oj). 
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The matrix of power spectra, S — (ui) J , is then obtained as 

g(u) =M _1 (cj)l(cj)M t " 1 (w). (C12) 
The diagonal elements of S are known as the power spectra, P a (ui). We find 

{fJ.%+w 2 )(a M f[x* P } +/J, M x* M ) + (a P x* M + npx*p)(a M f[x*p}\K(uj)\) 2 



(MmMp - - a M a P f'[x P ]Re[K (cj)]) 2 + (lu(^ m + (ip) - a M a P f [x* P ]lm[K (uj)]) 2 
a 2 p(a M f[xp} + Hmx" m ) + (fi 2 M + oj 2 )(a P x* M + fJ,px* P )) 



(C13) 



P P (u) = 

(mm^p - - aA./ap/'[a;p]Re[.?r(w)]) 2 + (u)(fj, M + Mp) - a M apf'[x P ]Im[K(uj)]) 2 
For constant delay, Jf(r) = S(r), we have ^(w) = e~ %uir , and we recover the result of [TU] . 

Appendix D: SIR-model with delayed recovery 

1. Gaussian approximation 

In the SIR-model with delayed recovery [22j . defined in the main paper, there are two independent types of 
particles, S and /. The number of recovered individuals follows from the constant overall population size. We 
write 

*=(?)• (Dl) 

In the model there are three reactions, two of which are delay reactions, labelled i — 2,3. The reaction rates 
are 

f f,(l-S-I)\ 

r[x(i)] = xPSI . (D2) 

V (i-x)0Si J 

The immediate effects of the reactions (i.e. at the time they are triggered) are given by 

1 0\ 

-11, (D3) 
1 1/ 

and the delayed effects on particle numbers are described by 

/0 \ 

m = I -1 I . (D4) 

The distributions of delay times are K 2 = K(t) and K 3 = Q(t), with K{t) and Q(t) as defined in the main 
paper. It turns out to be convenient to define K(t) — \K(t) and Q(t) = (1 — x)Q( T )- 
Putting all this together, and applying Eqs. (3, 4, 5) of the main paper we have 



F s (t, x) = -0S(t)I(t) + Ml ~ S(t) - J(t)) + p j dr Q{r)S{t - r)I(t - r), 

F T {t, x) = pS(t)I(t) -P [ dr [Q(t) + K( T )]S{t - r)I(t - r), (D5) 
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as well as 



B s ,s(t, t', x) = | \pS(t)J(t) + m(1 - S(t) - I(t)) +13 J dr Q(r)S(t - r)I(t - r)J S(t - t') 

-0Q{t - t')S(t')I{t') - /3Q(t' - t)S(t)J(t) J, 
B S!l {t,t',x) = U-pS(t)I(t)- p j dTQ(T)S(t-T)I(t-r)\5(t-t') 

+PQ{t - t')S(t')I(t') + p(p(t' -t) + K(t' - tf)S(t)I(t)\, 

B I>S {t,t',X) = U-PS(t)I(t)- p J dTQ( T )S(t-T)I(t-T)\S(t-t') 

+0Q{t' - t)S(t)I(t) + p(Q(t - t') + K{t - i')) S{t')I(t')y 
B ItI {t,l/,x) = | (pS(t)I(t)+p J dr (Q{T)+K(r)^jS(t-T)I(t-T)\s(t-l/) 

-p(Q(t - t')+K(t - t')^S(t r )I(t') ~ p(Q(t' - t) + K(t' - tj)S(t)I(t)j. (D6) 

2. Linear- noise approximation and spectra of quasi-cycles 

Concentrating again on parameter ranges in which the deterministic model, obtained for N — ¥ oo, has a fixed 
point (S* , I*), we decompose S = S* +£,s/VN and / = I* +£i/VN. Wi thin a sys tema tic expansion in powers 



of TV-i/2 t his allows us to substitute S(t) -> S* and I(t) -> I* in Eqs. (D5| and (D6) 



As in the model of gene regulation we next carry out a Fourier transform, and find 

MM€M = 3( W ) (D7) 

where 

MM = iwl-5(w,x*). (D8) 

For the SIR-model we have 

M(uj) = ( lUJ + " + J [1 ~ Q^ 1 * * + ft 1 ' Q W \ mq) 

— ' \ -P[l - K{u) - Q{u)]I* iw - P[l - K(u) - Q(u)]S* J ' V ' 

The power spectrum for the infectives, /, is given by Pi{oS) = ^|£/(cj)| 2 ^, and is found as 

Pi(") = |dct m^w (\Mss(oj)\ 2 B n (uj) - M ss (u)Mf s (u;)B IS (uj) 

-M is (uj)M* ss (lj)Bsi(lo) + \M si (uj)\ 2 B ss (lu)). (D10) 

For a specific instance of the model, i.e. for a specific choice of the delay kernel H (r) it is then a matter of first 
finding K(r), Q(t) and x> see main paper. A Fourier transform then gives K(u)), Q(u). The deterministic fixed 
point, S* and /*, can be found analytically for many kernels, as can the Fourier transforms, K{ui) and Q(uS) . 



Eq. (D10) can then be evaluated to obtain an analytical prediction for the power spectrum of quasi-cycles of 



the number of infectives in the population within the linear-noise approximation. 
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